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Abstract 


This report is the documentation of the work performed for the Subsonic Rotary Wing Project 
under the NASA’s Fundamental Aeronautics Program. It was funded through Task Number 
NNC10E420T under GESS-2 Contract NNC06BA07B in the period of 10/1/2010 to 8/31/2011. 
The objective of the task is to provide support for the development of variable speed power tur- 
bine technology through application of computational fluid dynamics analyses. This includes 
work elements in mesh generation, multistage URANS simulations, and post-processing of the 
simulation results for comparison with the experimental data. The unsteady CFD calculations 
were performed with the TURBO code running in multistage single passage (phase lag) mode. 
Meshes for the blade rows were generated with the NASA developed TCGRID code. The CFD 
performance is assessed and improvements are recommended for future research in this area. For 
that, the United Technologies Research Center's 1 V 2 stage Large Scale Rotating Turbine was se- 
lected to be the candidate engine configuration for this computational effort because of the com- 
pleteness and availability of the data. 


1 Contact information: waiming.to@nasa.gov (216-433-5937) 
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Introduction 


Unsteady multistage CFD code is a common tool for analyzing turbomachinery flows such as 
those encountered in the variable speed power turbine. Among the many choices of such flow 
codes, TURBO is chosen to be evaluated in this work. Since the code is largely employed in the 
analysis of compression systems, it is desirable to first focus on turbine aerodynamics before at- 
tempting for the more difficult heat transfer problem. Its performance in turbine application will 
be assessed and possible improvements are identified for future CFD simulation and TURBO 
code development. Also considered in this report is to provide some kind of documentation for 
best practices in performing turbine CFD simulations, particularly the procedure for starting 
from scratch with TURBO. 

A relevant multistage turbine configuration for the computation is the UTRC Large Scale Rotat- 
ing Rig (LSRR). Detail geometry information and aerodynamic data are available publicly from 
the UTRC reports [1,2]. A schematics of the rig copied from [1] is reproduced here in Figure 1. It 
is configured as a IV 2 stage turbine with blade counts of 22-28-28. Different configurations were 
reported and tested in this rig, but the particular configuration for the current CFD investigation 
is the IV 2 stage setup with the 50% nominal Statorl/Rotor gap because of the availability of the 
aero data from the experiment. The rotor tip has a clearance of 1.5% of the span. The flow pa- 
rameters chosen for the CFD is the 0.78 flow coefficient and 0.5% inlet turbulence intensity. The 
inlet Mach number is approximately 0.07 at this flow coefficient with a shaft speed of 410 rpm. 

Coordinates for each blade were published in the reports [1,2] for 3 spanwise positions (2%, 

50%, and 98%). This information is used to generate computational meshes for the blade pas- 
sages using the mesh generation code TCGRID [3]. H-mesh topology is employed through out 
the work. The computational domain begins at about 1 axial cord length of Stator 1 from its lead- 
ing edge and terminates at approximately 1 axial cord length of Stator 2 from its trailing edge 
(see Figure 2). 

The computations were performed on the NASA Glenn Research Center HX Compute Cluster. 
Each passage is divided into 16 blocks. A total of 48 processors were requested for each run in 
these simulations with each block assigned to a processor. 

Isolated Stator 1 calculations were also performed mostly for gauging the adequacy of grid qual- 
ity and numerical input parameters. Its computational exit boundary is extended to about 1 axial 
cord downstream instead of ending at the Stator 1 /Rotor interface. 
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In the next section, details on generating the meshes and preprocessing of the data for TURBO 
are described. The Computation Details section provides additional information about setting up 
the calculations, computing requirements and procedures of the simulations. Following that, the 
results from these calculations are compared to the test data reported in [2]. Finally, concluding 
remarks are given in the last section. 

An accompanying DVD is also available upon request an approval by NASA. It contains all the 
files necessary for restarting the calculations for the isolated Stator 1 and the V /2 stage. In addi- 
tion, the TURBO source code, the airfoil coordinates and their meshes are also available from the 
same DVD. 


Mesh Generation and Preprocessing 


Blade profiles for each of the blade are given in [1,2] as coordinates at 3 different radial posi- 
tions. They are digitally reproduced here and plotted in Figure 3 to Figure 5. This data are pro- 
vided in the DVD under the AirfoilCoordinates folder. Since the resolution around the leading 
and trailing edges from the raw data is insufficient for making acceptable meshes, these regions 
are fitted with circles to produce a smoother definition of the blade contour. Radii of these circles 
and coordinates of their centers are shown in the figures. The blade contours are given in inches 
at the 2%, 50%, and 98% span corresponding to the hub, midspan, and tip profile denoted in the 
figures. 

This information is utilized in the input to TCGRID for generating the computational mesh for 
each blade passage. The radial grid lines at the blade row interface have to be matched up from 
both sides. It is an implicit assumption made in the TURBO code for interpolation across the in- 
terface. Input file format and input parameters for running TCGRID are well documented in [3]. 
The input files for generating the meshes are given in the DVD under the Mesh folder as 
tcgrid. SI. input, tcgrid.Rl. input, and tcgrid.S2. input for Statorl, Rotorl, and Stator2 respectively. 
H-mesh topology is employed for the blade passages. The rotor mesh has 5 radial points in the 
tip clearance region. The mesh output from TCGRID is written in PLOT3D single grid format in 
Little Endian which is then converted to the multi grid Big Endian format in order to be read in 
by GUMBO for further mesh preprocessing. Care on the distance of the first grid point from the 
wall should be exercised as the turbulence model used in TURBO is of a high Reynolds number 
form of d-s model with wall damping function. The assessment of y + is done with a series of pre- 
liminary runs for an isolated Stator 1 with different near wall mesh density. While the range of y + 
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for the turbulence model to be valid lies between 20 and 100, a limited number of points with y + 
< 20 in the flow is permitted in the computation. The final revision of the meshes for the simula- 
tions are shown in Figure 2, with mesh densities of (206, 64, 61), (216, 64, 61), and (216, 64, 61) 
for Stator 1 , Rotor, and Stator 2 respectively. 

Mesh and block manipulations, block boundary condition and block connectivity settings, and 
blade row interface information are conveniently set up and controlled with GUMBO, a graphi- 
cal user interface tool for manipulating the initial data and grid. Initial flow field, if available can 
also be set up with GUMBO simultaneously. However, this is only restricted to APNASA data 
type format. Since there is no initial flow field available, the computation will have to start from 
scratch with uniform flow. Files for the block connectivity, block boundary conditions, and blade 
row interfaces, written by GUMBO are all that is required for starting a TURBO computation. It 
should be remarked that these files can be generated manually without necessarily going through 
GUMBO. In addition, a set of input files are also needed by TURBO for setting up parameters 
such as inlet and exit boundary conditions, reference conditions, shaft speed, time step size and 
other numerical parameters. Details of the description of input parameters can be found from the 
link: http://www.simcenter.msstate.edu/docs/msu_turbo/ . The list of files created from this pre- 
processing step is provided under the CASE directory in the DVD. 

Initially each blade row passage is divided into 5 blocks. Although this block density resulted in 
a reasonable turnaround time for the isolated Stator 1 simulation, it becomes evident that larger 
block density is needed for stage computation in order to speed up the calculations because of the 
low Mach number flow. The unsteady nature of the blade row interaction has created a long nu- 
merical transient despite of a peak local Mach number of about 0.26 on the suction side of the 
Stator 1. Eventually, a 16 blocks/passage is adopted for the subsequent simulations. Each block 
is assigned to a processor and will need about 300MB of memory which is well within HX Clus- 
ters’ 1GB per processor limit. For the stage computation, a total of 48 processors are required. 


Computation Details 


The calculations were carried out with 60 time steps per rotor passing. The solution normally 
converges in 10 to 15 revs which is equivalent to a wallclock time of 24 to 36 hours on the HX 
Clusters. For computation starting from scratch, additional time is needed for the flow field to 
arrive at a reasonable level. This is because of the necessity for very small time step sizes at the 
early stages of the calculation. 
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For the 22-28-28 blade count, performing a periodic sector computation of 11-14-14 blade count 
(half wheel) would have required 624 processors while maintaining the same computational effi- 
ciency of 16 blocks/passage. That is beyond the current HX limit. The calculations will have to 
be performed on a larger computer such as the NAS compute clusters which has other constraints 
to consider. The alternative is to reduce the number of blocks per passage by half by sacrificing 
the wallclock time, resulting in roughly doubling the turnaround time. Based on these considera- 
tions, it is decided to run TURBO on HX in single passage phase lag mode instead. 

Adiabatic wall condition is assumed so there will be no heat transfer calculation. Uniform total 
condition with zero flow angles is imposed at the inflow boundary. This assumes no boundary 
layer present from the inlet. At the outflow boundary, radial equilibrium flow is assumed with 
static pressure (back pressure) imposed at the hub. Its level is to be adjusted during calculations 
in such a way until a desired flow coefficient is reached. It is a tedious and time consuming proc- 
ess involving restarting TURBO and sometimes resetting the time information of initial flow 
field. After running a number of trial cases, the exit static pressure is found to be somewhere be- 
tween 95% and 96% of the inlet total pressure for the VA stage. 

As mentioned earlier, the inlet turbulence intensity for this case was measured at 0.5% in the ex- 
periment. However, the turbulence kinetic energy, &, and the dissipation rate, s, are required as 
the boundary conditions at the inflow boundary. While L can be reasonably approximated by 1.5 
times the turbulence intensity assuming turbulence isotropy, s is to be estimated by assuming an 
eddy viscosity at the inflow boundary. Trial runs with eddy viscosity values of 5, 50, and 300 
show insignificant sensitivity on the midspan pressure coefficients of the blade. As a result, a 
value of 50 is used for the present CFD analyses. 

Besides back pressure, there are instances when the time step size or the shaft speed needs to be 
modified in the middle of the calculation, and the same procedures for restarting the calculation 
will be followed. This process is not as bad as starting from scratch since the restart files or the 
flow solution files from some previous TURBO runs, converged or not, are usually available. 
Nevertheless, it will still have to undergo a numerical transient period. 

Starting the computation from scratch is tedious but rather straightforward, the flow is initialized 
with either an uniform axial flow or throughflow. It is necessary to begin with small initial time 
step size or the CFL number and monitor the progress by gradually dialing up the time step size 
or the CFL number to the desired level. In some instances with high speed flow, it may even be 
necessary to lower the shaft speed in order to stabilize the numerical transient. For the current 
case, 500 steps per rotor passing is used initially. 
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A final remark on speeding up the computation by raising the inlet Mach number: it is found that 
the rate of convergence was not significantly improved with an inlet Mach number of about 0.15. 
In this case, the shaft speed has to be raised to about 900 rpm in order to maintain a similar flow 
coefficient. This resulted in a peak local Mach number of about 0.6 (compared to 0.26 with Moo = 
0.07) on the suction side of Stator 1, which can still be considered as having negligible effect on 
the blade midspan pressure distribution. The exercise was carried out in hope of improving the 
overall wallclock time by elevating the inlet Mach number. It turns out that the combination of 
the numerical and computational parameters such as Newton iteration number, Jacobian matrix 
update frequency, Gauss-Seidel iteration number, time step size, and number of blocks per pas- 
sage have a more significant impact on the computational speed, particularly the latter one. 


Results 


Most of the computations reported in this work converged in 10 revolutions or less based on the 
time averaged inlet mass flow. Once converged, time averaged flow field is obtained by averag- 
ing the time record of the unsteady field over the last 2 revs. The flow field consisting of the 5 
dependent variables (q, qu, qv, qw, e) are then used to compute other state variables or flow 
quantities such as pressure or velocity through either the equation of state or the isentropic rela- 
tions. The assumption made in this approach is essentially that the quantity derived from the time 
averaged dependent variables is the “same” as the time averaged of the unsteady derived quan- 
tity. The results presented here are all manipulated in this way. 

Constant specific heat ratio is assumed because of the small temperature change in the flow field. 
Inlet total pressure and temperature are maintained at 101,325 Pascal and 294°K. Static pressure 
at the hub of the exit plane is adjusted until the desired flow coefficient is reached. The rotor tip 
clearance gap region is not gridded, but rather a clearance model is used in the calculations. In 
the subsequent paragraphs, results from the isolated Stator 1 calculation are presented first along 
with some observations when comparing to the stage Stator 1 results, and then followed by the 
discussion of the stage results. 

A typical mass flow time history is shown in Figure 6. The time averaged inlet mass flow for 
both the isolated Stator 1 and the 1 Vi stage CFD results are plotted on the same figure for com- 
parison. The isolated Stator 1 computation is performed with a similar mesh distribution used for 
Stator 1 in the stage computation, except for extending the exit plane to about 1 axial cord length 
downstream of the trailing edge. The exit static pressure is adjusted in a similar way as described 
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earlier for the stage until the resulting flow coefficient is about 0.78 which is the target flow coef- 
ficient for the stage calculation. It should be mentioned that the isolated Stator 1 adjustment is 
less challenging because of its steady nature and faster convergence. As seen from the figure, the 
stage result displays a higher level of fluctuation probably due to the unsteady nature of blade 
row interactions. The isolated Stator 1 result shows that the calculation converges very quickly 
and that the flow is steady compared to the 1 Vi stage result. 

Figure 7 is a plot of the time averaged pressure coefficient for the isolated Stator 1. The circles 
represent test data from [2] for the VA stage experiment and are shown here for reference. The 
pressure coefficient, C p , is defined as (Po - pjMQoVRm 2 , where Po and qo are the averaged inlet 
total pressure and density at midspan respectively, p is the blade surface static pressure at mid- 
span, and VRm is the reference speed which is assumed to be the midspan rotational speed of the 
rotor in the stage configuration. The data has an uncertainty in measurements of ±2% for C p and 
±1% for the flow coefficient. The result shows that the pressure coefficient is not sensitive to 
whether Stator 1 is configured as an isolated row or in a stage (see Figure 10a for Stator 1 pres- 
sure coefficient from the stage calculation). Results from both stator only and stage calculations 
show that the blade surface pressure and the wake in the midspan region is two dimensional in 
nature. 

Figure 8 shows the time averaged total pressure contour for the isolated Stator 1 and the VA stage 
at midspan. The scale of the contour for each of the blade passages is shown next to it. The grey 
vertical line shown behind the blade for the stage contour in Figure 8(b) is where the circumfer- 
ential survey were made in the experiment along the midspan. No flow separation is observed 
from the calculations. It is seen from the figure that Stator 1 wake for the stage result is thicker 
than that for the isolated Stator 1 . Upon inspection of the mesh, it is probably due to the large 
numerical diffusion from the poor mesh quality in the inter-blade row gap that contributes to the 
wake blooming. Comparing the 2 Stator 1 meshes, the isolated Stator 1 mesh changes more 
gradually extending to the exit plane over a longer distance, whereas the stage Stator 1 mesh 
must terminate at the blade row interface resulting in large stretching and skewing. 

Pressure coefficients and blade exit profiles for the stage CFD results are shown in Figure 9 to 
Figure 12 together with the test data from [2]. The computational results are time averaged over 
the last 2 revs. Midspan pressure coefficient for Stator 1, Rotor 1, and Stator 2 are shown in Fig- 
ure 9a, 9b, and 9c respectively. The pressure coefficient, C p , is defined in the same way as be- 
fore. The circles from the figures represent the test data. CFD results from 2 different flow coef- 
ficients are shown in the same plot. The flow coefficients are not exactly at the experimental 
condition of 0.78 primarily because of the mass flow sensitivity to the back pressure adjustment 
while converging the unsteady flow to the target time averaged flow coefficient. Nevertheless, 
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one can deduce from the trend in Figure 9 that higher flow coefficient results in higher pressure 
coefficient. For Stator 1, that would have made the CFD results closer to the test data had a 
higher flow coefficient been achieved. However, that would make the pressure coefficient drift- 
ing away from the data for Rotor 1 and Stator 2, making it worst on the suction side in compari- 
son. 

Another comment from Figure 9 is that the CFD result is progressively worst for the subsequent 
blade rows, particularly on the suction side of the blade although the shape of the pressure coef- 
ficient profile follows the test data reasonably close. It is not clear if such discrepancy is attrib- 
uted to issues related to the Rotor 1 and Stator 2 blade geometry and their meshes, or the poor 
Stator 1 wake upon the rotor. 

Circumferential variations in flow speed, flow angle, and static pressure at midspan behind Sta- 
tor 1 are shown in Figure 10. The results are presented at a distance of 17% axial cord from the 
its trailing edge. The solid lines represent the CFD results in the figures. The test data are shown 
in dots which are copied from the figures reported in [2]. Similar plots for Rotor 1 and Stator 2 
are illustrated in Figure 1 1 and Figure 12 respectively, except that the rotor data are presented at 
a distance of 36% axial cord from its trailing edge and the stator at a distance of 14% axial cord 
from its trailing edge, shown as grey lines in Figure 8 where experimental measurements were 
made. Moreover, the flow speed and flow angle are plotted in absolute frame of reference for the 
stators and in relative frame for the rotor. The static pressures are normalized in the same way as 
the pressure coefficients and the flow speed is normalized by VR m as before. 

The comparisons show that only qualitative agreement with the test data can be drawn from 
these CFD results. For example, all 3 exit profiles for the flow speed in Figure 10 to Figure 12 
indicate that the wakes from their respective blades are too diffused. Nonetheless, their average 
values are generally not too far off from the measurements except for the rotor exit flow angle 
which is under- turned by 6°. The following table shows the time and circumferentially averaged 
exit flow angles comparison with the data. Note that the angles given in the table are measured 
form the vertical, whereas those in the figures are measured from the horizontal, following the 
same convention from the UTRC report [2]. 

A closer examination reveals that the worst mesh quality is usually found in the inter-blade row 
gap regions or near the leading and trailing edges. The H-mesh topology for large turning angle 
turbine blade frequently results in large mesh expansion ratio and sheared cells where numerical 
diffusion are substantial. 
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Axial Position 


Test Data 


CFD 


1 7% cord from Stator 1 TE 

67.6° (Absolute) 

66.33° (Absolute) 

36% cord from Rotor 1 TE 

58.3° (Relative) 

52.30° (Relative) 

14% cord from Stator 2 TE 

66.3° (Absolute) 

65.82° (Absolute) 


Table 1 . Time and Circumferentially Averaged Flow Angles Downstream of Each Blade Row. 
(Axial Position Given in Percent of Axial Cord From the Trailing Edge) 


Conclusion 


Blade contours for the UTRC’s IV 2 stage LSRR in turbine configuration were generated from the 
coordinates given in [1,2]. The leading and trailing edges were fitted with circles to generate bet- 
ter definition blade profiles. This information was then used as input to TCGRID in generating 
meshes for the blade passages. The unsteady calculations were performed with TURBO in single 
passage phase lag mode. Time averaged midspan pressure coefficient for each blade are com- 
pared with the test data. Time averaged flow speed, flow angle, and static pressure at midspan 
blade exit plane as described in [2] are computed and compared with the data. 

The comparisons show that the CFD results are only in qualitative agreement with the test data 
and indicate that the poor mesh quality might be the cause for the discrepancies. This is an inher- 
ent drawback for turbine blade gridded with H-mesh. Ultimately, a mesh topology such as an O- 
H mixed type of mesh or a multi-zone type is more appropriate and likely to be more successful 
for turbine blade CFD simulations. 

For this to happen, supports in mesh generation are needed for sophisticated multistage mesh to- 
pology. This necessitates further TURBO code development and support in the area of general 
block to block MPI communication. Although built for that in mind, TURBO has not been rigor- 
ously tested in this aspect. Therefore, additional efforts in testing and validation in supporting 
this feature is recommended. 
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Figures 



Figure 1. UTRC Large Scale Rotating Rig Schematics (diagram taken from [1] for illustration). 
The Present CFD Analysis is Configured for the 50% Nominal Sl/Rl Axial Gap. 
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Figure 2. H-Mesh for the VA Stage LSRR Computation. 
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LSRR Stator 1 Hub Profile LSRR Stator 1 Midspan Profile 

Pitch 6.88865 inches Pitch 7.7 1 1 1 8 inches 




LSRR Stator 1 Tip Profile 


Pitch 8.53371 inches 



Figure 3. Stator 1 Blade Contours. 
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LSRR Rotor 1 Hub Profile 


LvSRR Rotor 1 Midspan Profile 


Pitch 5.41251 inches 


Pitch 6.05879 inches 




LSRR Rotor 1 Tip Profile 


Pitch 6.70506 inches 



Figure 4. Rotor 1 Blade Contours. 
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LSRR Stator 2 Hub Profile LSRR Stator 2 Midspan Profile 

Pitch 5.41251 inches Pitch 6.05879 inches 




LSRR Stator 2 Tip Profile 


Pitch 6.70506 inches 



Figure 5. Stator 2 Blade Contours. 
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Figure 6. Mass Flow Convergence History Comparison. 
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Figure 7. Isolated Stator 1 Midspan Pressure Coefficient. 
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(a) Isolated Stator 1 Results 



(b) IV 2 stage Results 


Figure 8. Midspan Time Averaged Total Pressure Contour. 
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LSRR Stator 1 Midspan Pressure Coefficeint 


Test Data at Flow Coefficient of 0.78 (+/- 1%) 



Figure 9a. Stator 1 Midspan Time Averaged Pressure Coefficient. 


LSRR Rotor 1 Midspan Pressure Coefficeint 


Test Data at Flow Coefficient of 0.78 (+/- 1%) 



Axial Cord Distance 


Figure 9b. Rotor 1 Midspan Time Averaged Pressure Coefficient. 
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LSRR Stator 2 Midspan Pressure Coefficeint 


Test Data at Flow Coefficient of 0.78 (+/- 1%) 



Figure 9c. Stator 2 Midspan Time Averaged Pressure Coefficient. 
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Figure 10. Midspan Stator 1 Exit Profile at 17% Axial Cord From Trailing Edge. 
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Figure 11. Midspan Rotor 1 Exit Profile at 36% Axial Cord From Trailing Edge. 
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Figure 12. Midspan Stator 2 Exit Profile at 14% Axial Cord From Trailing Edge. 
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